Data Analysis and Visualization
Paper 2: Community Data for 4 countries
WORK IN PROGRESS

Author

Gagan Atreya

Published

November 21, 2023

Section 1. Data cleaning and Demographics

Display code
options(digits = 2)
rm(list=ls())

## Install "pacman" package if not installed
# (remove the # symbol from the line below):
# install.packages("pacman")

## Load R packages:
pacman::p_load(data.table, tidyverse, haven, labelled, vtable, 
               psych, scales, weights, clipr, forcats,
               stargazer, ggthemes, ggcharts, geomtextpath,
               corrplot, tm, readxl, patchwork, modelsummary,
               gt, lme4, car, lmerTest, 
               ggeffects, magrittr, broom, broom.mixed,
               backports, effects, interactions, plyr, sjPlot)

## Import datasets:
## These are newer datasets with new variables created in the individual analyses:
dsgmb <- fread("~/Dropbox/2023_Gagan/Paper2/data/dsgmb.csv")
dsgmb$Country <- "Gambia"
dsgmb$id <- 1:nrow(dsgmb)
dsgmb$ID <- paste0("GMB",dsgmb$id)

dspak <- fread("~/Dropbox/2023_Gagan/Paper2/data/dspak.csv")
dspak$Country <- "Pakistan"
dspak$id <- 1:nrow(dspak)
dspak$ID <- paste0("PAK",dspak$id)

dstza <- fread("~/Dropbox/2023_Gagan/Paper2/data/dstza.csv")
dstza$Country <- "Tanzania"
dstza$id <- 1:nrow(dstza)
dstza$ID <- paste0("TZA",dstza$id)

dsuga <- fread("~/Dropbox/2023_Gagan/Paper2/data/dsuga.csv")
dsuga$Country <- "Uganda"
dsuga$id <- 1:nrow(dsuga)
dsuga$ID <- paste0("UGA",dsuga$id)

## Correct asterisk pattern for stargazer tables:
starpattern <- "<em>&#42;p&lt;0.05;&nbsp;&#42;&#42;p&lt;0.01;&nbsp;&#42;&#42;&#42;p&lt;0.001</em>"

## List of variables to retain from all datasets:
list1 <- c("ID", "Country", "age", "gender", 
           "ses", "jobnature", "religion", "married", "education", "Ethnicity",
           "IGF01", "IGF02", "IGF03", "IGI01", "IGI02", "IGI03", 
           "OGF01", "OGF02", "OGF03", "OGI01", "OGI02", "OGI03",
           "ENDBCL01", "ENDBCL02", "ENDBCL03", "ENDBBL01", "ENDBBL02", "ENDBBL03",
           "EXPBCL01", "EXPBCL02", "EXPBCL03", "EXPBBL01", "EXPBBL02", "EXPBBL03", 
           "empathic_concern_01", "empathic_concern_02", "empathic_concern_03",
           "perspective_taking_01", "perspective_taking_02", 
           "perspective_taking_03", "perspective_taking_04", "history_discrimination",
           "og_hostility", "og_cooperation", "fight_outgroup", "imagistic",
           "event_positive_affect", "event_negative_affect", "event_episodic_recall",
           "event_shared_perception", "event_event_reflection", 
           "event_transformative_indiv", "event_transformative_group",
           "freq_positive_contact", "freq_negative_contact", 
           "sprf", "exp_religious_freedom")

## Subset datasets to only the columns in the above list:
dsgmb1 <- dsgmb[, ..list1]
dspak1 <- dspak[, ..list1]
dstza1 <- dstza[, ..list1]
dsuga1 <- dsuga[, ..list1]

## Merged dataset with needed columns only
ds <- rbind(dsgmb1, dspak1, dstza1, dsuga1)

## Rename the "Event_" columns with title case:
ds01 <- ds[, 1:44]
ds02 <- ds[, !(2:44)]
colnames(ds02) <- stringr::str_to_title(colnames(ds02))
ds02$ID <- ds02$Id
ds02 <- ds02[, !1]
ds <- merge(ds01, ds02, by = "ID")
rm(ds01, ds02)

Data cleaning

Criteria:

  • Survey duration:

not available for this dataset as it was manual interviews

  • Bad open ended responses:

not applicable since all are manual responses

  • SD on unrelated scales

  • INCOME VARIABLE

some cases with outlier responses in “income” variable need to be removed. THIS PART NEEDS TO BE REDONE

  • Missing data in scales

Missing data:

Display code
mds <- map(ds, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       2       2       3       2      34 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data")+
  theme_bw()

lp02

Display code
## Individual country datasets:
dsgmb <- ds[ds$Country == "Gambia", ]
dspak <- ds[ds$Country == "Pakistan", ]
dstza <- ds[ds$Country == "Tanzania", ]
dsuga <- ds[ds$Country == "Uganda", ]

Missing data: Gambia

Display code
mds <- map(dsgmb, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.4     0.8     1.3     1.7     9.3 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Gambia")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Gambia")+
  theme_bw()

lp02

Missing data: Pakistan

Display code
mds <- map(dspak, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       2       0     100 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp02

Pak after removing BCL01:

Display code
dspak22 <- subset(dspak, select = -c(EXPBCL01))

mds <- map(dspak22, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00    0.00    0.07    0.00    0.79 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:54,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Pakistan")+
  theme_bw()

lp02

Missing data: Tanzania

Display code
mds <- map(dstza, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       5       6       7       8      44 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Tanzania")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Tanzania")+
  theme_bw()

lp02

Missing data: Uganda

Display code
mds <- map(dsuga, ~mean(is.na(.))*100) 
mds <- as.data.frame(mds)
mds1 <- as_tibble(cbind(variable = names(mds), t(mds)))
mds1$percent_missing <- mds1$V2

mds1$percent_missing <- as.numeric(mds1$percent_missing)
summary(mds1$percent_missing)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.7     0.9     1.1     1.3     5.1 
Display code
mds2 <- select(mds1,-c(2))
mds2 <- as.data.frame(mds2)
mds2$percent_missing <- as.numeric(mds2$percent_missing)
mds2 <- as.data.table(mds2)
mds2a <- mds2[1:30,]
mds2b <- mds2[31:55,]

lp01 <- mds2a %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Uganda")+
  theme_bw()

lp01

Display code
lp02 <- mds2b %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = variable,
                 y = percent_missing,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "% missing data: Uganda")+
  theme_bw()

lp02

Demographic Variables:

Display code
ds$gender <- ifelse(ds$gender == "Male", "Male",
             ifelse(ds$gender == "Female", "Female", NA))

## ses:
ds$ses <- ifelse(ds$ses == "", NA, ds$ses)

## jobnature:
ds$jobnature <- ifelse(ds$jobnature == "", NA, ds$jobnature)

#sentence case for jobnature:
ds$jobnature <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$jobnature, perl = TRUE)

ds$jobnature <- ifelse(ds$jobnature == "Non-government/self-employed", 
                       "Non-government", ds$jobnature)

## Religion:
ds$religion <- ifelse(ds$religion == "", NA, ds$religion)

ds$religion <- ifelse(ds$religion == "Christian (Catholic)", "Christian: Catholic", 
               ifelse(ds$religion == "Christian (Protestant)", "Christian: Protestant",
               ifelse(ds$religion == "Muslim (Shia)", "Muslim: Shia",
               ifelse(ds$religion == "Muslim (Sunni)", "Muslim: Sunni", ds$religion))))

## Marital status

ds$married <- ifelse(ds$married == "Not married", "Unmarried", ds$married)
ds$married <- ifelse(ds$married == "", NA, ds$married)

## Individual country datasets:
dsgmb <- ds[ds$Country == "Gambia", ]
dspak <- ds[ds$Country == "Pakistan", ]
dstza <- ds[ds$Country == "Tanzania", ]
dsuga <- ds[ds$Country == "Uganda", ]

Demographic variables in the analysis:

  • Age
  • Gender
  • Socio-economic status
  • Nature of employment
  • Religious affiliation
  • Marital status
  • Ethnicity

Variable: Sample size by Country

Display code
tbl01 <- table(ds$Country)

## Table of user language by country:
tbl01

  Gambia Pakistan Tanzania   Uganda 
     237      505      352      448 
Display code
## Sample size by country:

lp01 <- ds %>% 
  # drop_na(Country) %>%
  lollipop_chart(x = Country,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Sample size by country")+
  theme_bw()

lp01

Variable: Age

Display code
summary(ds$age)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
     18      26      33      37      45      92      29 
Display code
ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 37.00", 
                 xintercept = 37.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Age", 
       y = "Frequency", 
       title = "Age distribution (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  labs(x = "Age", 
       y = "Frequency", 
       title = "Age distribution by country")+
  facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Gender

Display code
ds$gender <- ifelse(ds$gender == "Male", "Male",
             ifelse(ds$gender == "Female", "Female", NA))

lp02 <- ds %>% 
  drop_na(gender, age, Country) %>%
lollipop_chart(x = gender,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Gender distribution (full sample)")+
#  facet_wrap(~Country, nrow = 2)+
  theme_bw()

lp02

Display code
bp01 <- ds %>% drop_na(gender, age) %>% 
  ggplot(aes(y = age, 
             x = gender))+
geom_boxplot(fill = "grey")+
  labs(y = "Age",
       x = "",
       title = "Age and gender distribution by country")+
  facet_wrap(~Country, nrow = 2)+
  coord_flip()+
  theme_bw()

bp01

Variable: Socio-economic status

Display code
ds$ses <- ifelse(ds$ses == "", NA, ds$ses)
table(ds$ses)

Lower middle       Middle        Upper Upper middle 
         479          831           27          172 
Display code
ds %>% drop_na(ses) %>%
lollipop_chart(x = ses,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Socioeconomic status (full sample)")+
  theme_bw()

Variable: Nature of employment

Display code
table(ds$jobnature)

    Government Non-government          Other        Retired  Self-employed 
           329            295              2              1            773 
       Student 
             6 
Display code
ds %>% drop_na(jobnature) %>%
lollipop_chart(x = jobnature,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Nature of employment (full sample)")+
  theme_bw()

Variable: Religious affiliation

Display code
table(ds$religion)

    Christian: Catholic        Christian: Other  Christian: Pentecostal 
                    460                      43                      33 
  Christian: Protestant                   Hindu           Muslim: Other 
                    209                      37                       4 
           Muslim: Shia           Muslim: Sunni                   Other 
                    188                     433                      87 
                   Sikh Spiritual not religious 
                      9                      14 
Display code
lp05 <- ds %>% drop_na(religion) %>%
lollipop_chart(x = religion,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Religious distribution (full sample)")+
  theme_bw()

lp05

Variable: Marital status

Display code
table(ds$married)

  Married     Other Unmarried 
      957        94       483 

Variable: Education

Display code
table(ds$education)

                                         
                                      44 
                                       4 
                                       8 
                                   4m  4 
                                       1 
                                       7 
                                       7 
                         Advance Diploma 
                                       1 
     Advance Diploma in Business Studies 
                                       1 
           Advance Diploma in Management 
                                       3 
              Advance Diploma in Nursing 
                                       1 
   Advance Diploma in Tourism Management 
                                       1 
                           Advance level 
                                       1 
                                   B.S.C 
                                       1 
                                Bachelor 
                                     104 
                         Bachelor Degree 
                                       2 
     Bachelor in Business Administration 
                                       1 
                       Bachelor's Degree 
                                       3 
  Bachelor's Degree in Political Science 
                                       1 
                     Bachelors (ongoing) 
                                       1 
                        Bachelors Degree 
                                       1 
                  Bachelors in Education 
                                       1 
                       Bachelors of Edu. 
                                       1 
                                     Bsc 
                                       1 
                                     BSc 
                                       3 
                                     BSC 
                                       5 
                  BSc (Arabic Education) 
                                       1 
                   BSc (Hons) Management 
                                       1 
                       BSc in Accounting 
                                       1 
                         BSc in Economic 
                                       1 
                        BSc in Economica 
                                       1 
                        BSc in Education 
                                       1 
                   BSc in Human Resource 
                                       1 
                       BSc in Management 
                                       1 
                BSc in Political Science 
                                       1 
                      BsSc of Accounting 
                                       1 
                             Certificate 
                                       3 
   Certificate in International Relation 
                                       1 
       Certificate in Management Studies 
                                       1 
                       Certificate Level 
                                       1 
    Certificate on Animal Health Product 
                                       1 
                            Certtificate 
                                       1 
                        Chartered Banker 
                                       1 
                                   Cheti 
                                       6 
                                    Chuo 
                                       6 
                             Chuo Ckikuu 
                                       1 
                              Chuo kikuu 
                                      11 
                              Chuo Kikuu 
                                       3 
                             Claiasa nne 
                                       1 
                                 College 
                                       9 
                           College (HTC) 
                                       1 
                           College (HWD) 
                                       1 
                     College Certificate 
                                       1 
                         currently in s4 
                                       1 
                                Darasa 7 
                                       8 
                                Darasa 8 
                                       1 
                             Darasa la 7 
                                       1 
                          Darasa la Saba 
                                       1 
                              Darsa la 7 
                                       2 
                              Darsa La 7 
                                       3 
                           Darsa la saba 
                                       1 
                                  Degree 
                                      15 
               Degree in Islamic Studies 
                                       1 
                          Degree-Shahada 
                                       1 
                           Didato channe 
                                       1 
                           Digri/Shahada 
                                       1 
                                 Diploma 
                                      22 
       Diploma 2 Business Administration 
                                       1 
           Diploma in Community Policing 
                                       1 
                          Diploma in ICT 
                                       2 
                   Diploma in Management 
                                       3 
                     Diploma in Theology 
                                       1 
           Diploma in Travel and Tourism 
                                       1 
                                Diplomat 
                                       1 
                            Elimu Msingi 
                                       1 
           Form 3 (Old education System) 
                                       2 
           Form 3 (Old Education System) 
                                       1 
                                  Form 4 
                                      15 
                                  Form 5 
                                       1 
           Form 5 (Old education System) 
                                       2 
           Form 5 (Old Education System) 
                                       1 
                                  Form 6 
                                       3 
                               Form Five 
                                       2 
                             Form Five 5 
                                       2 
                               Form Four 
                                      73 
                             Form Four 4 
                                       3 
                             Form Four 5 
                                       1 
                                 Form IV 
                                       1 
                              Form One 1 
                                       1 
                                Form Six 
                                      14 
                              Form Three 
                                       1 
                                Form Two 
                                       6 
                              Form Two 2 
                                       1 
                    Gambia College (PTC) 
                                       1 
                       GCE Level(WASSCE) 
                                       1 
                  Grade 8 ( Upper Basic) 
                                       1 
                        Graduate Diploma 
                                       1 
                             High School 
                                       4 
                    High School (WASSCE) 
                                       1 
                 High School Certificate 
                                       2 
                       High School Level 
                                       1 
           High School-Islamic Education 
                                       1 
                          Higher Diploma 
                                       1 
       Higher Diploma in Education (HDC) 
                                       1 
        Higher Teacher Certificate (HTC) 
                                       7 
      Higher Teacher's Certificate (HTC) 
                                       5 
                            Hotel School 
                                       1 
                                     HTC 
                                       2 
                             HTC Primary 
                                       1 
                                 HTC/HNC 
                                       1 
                              Illiterate 
                                      10 
            Informal Education (Madrasa) 
                                       1 
                            Intermediate 
                                      85 
                       Isachetovs Degree 
                                       1 
                          Islamic School 
                                       1 
              Islamic Traditional School 
                                       1 
                                      IV 
                                       1 
                       Junior Sec School 
                                       1 
                        Junior Secondary 
                                       1 
                                Kidato 4 
                                       3 
                            Kidato cha 1 
                                       1 
                            Kidato cha 4 
                                       1 
                            Kidato cha 6 
                                       2 
                          Kidato cha nne 
                                       2 
                         Kidato Cha seta 
                                       1 
                           Kidato Channe 
                                       1 
                           Kidato chaune 
                                       1 
                                    La 4 
                                       3 
                                    La 7 
                                       7 
                                    LA 7 
                                       8 
                                  La nne 
                                       4 
                              La nne (4) 
                                       1 
                          La nne mkoloni 
                                       1 
                                 La saba 
                                      17 
                               La saba 7 
                                       8 
                           Lanne Mkoloni 
                                       2 
                                  Lasaba 
                                       1 
                                     LLB 
                                       1 
                                  Master 
                                     161 
                          Masters Degree 
                                       2 
                    Masters in Economics 
                                       1 
                    Masters in Sociology 
                                       1 
                             Matriculate 
                                      34 
Memorized the Quran (Informal Education) 
                                       1 
                                  Middle 
                                      53 
                                   MPhil 
                                      56 
                                     MSC 
                                       1 
                                  Msingi 
                                       3 
                                    None 
                                      50 
                                 O Level 
                                       3 
                         Ordinary Levels 
                                       1 
                                     PhD 
                                       3 
                           Post Graduate 
                                       2 
                   Post Graduate Diploma 
                                       1 
                                 Primary 
                                     213 
                          Primary School 
                                       3 
       Primary Teacher Certificate (PTC) 
                                       2 
     Primary Teacher's Certificate (PTC) 
                                       1 
                                     PTC 
                                       1 
            Quranic High School Graduate 
                                       1 
                          Quranic School 
                                       1 
                              Seconary 4 
                                       1 
                               Secondary 
                                     143 
                         Secondary  Four 
                                       1 
                         Secondary level 
                                       1 
                        Secondary School 
                                       6 
                                Sekonorw 
                                       1 
                Senior School Graduation 
                                       1 
                        Senior Secondary 
                                       2 
                 Senior Secondary School 
                                      16 
    Senior Secondary(WASSCE Certificate) 
                                       1 
                                 Shadada 
                                       3 
                                 Shahada 
                                       5 
                             Shahada yai 
                                       1 
                             Shule Msiry 
                                       1 
                         Shule ya msingi 
                                       1 
                                Sokondal 
                                       1 
                              StaShahada 
                                       2 
                                 Std VII 
                                       3 
                                 STD VII 
                                       8 
                                Tertiary 
                                      33 
                      Tertiary Education 
                                      10 
                          Tertiary level 
                                       1 
             Traditional Quranic Schools 
                                       1 
                            Udergraduate 
                                       1 
  Undergraduate Degree in Communications 
                                       1 
                              University 
                                      15 
                       University Degree 
                                       7 
                University of the Gambia 
                                       1 
                                 Uzamili 
                                       1 
                                      VI 
                                       1 
                                 VII Std 
                                       1 
                              Vocational 
                                       8 
                                  WASSCE 
                                      16 
          WASSCE examination certificate 
                                       1 

Variable: Ethnicity

Display code
## Gambia:
eth <- as.data.frame(table(dsgmb$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 2, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dsgmb2 <- dsgmb[, c("Ethnicity")]

dsgmb2$ethnicity <- ifelse(dsgmb2$Ethnicity %in% l1, dsgmb2$Ethnicity, "Other")
dsgmb2$ethnicity <- ifelse(dsgmb2$ethnicity == "Serere", "Serer",
                 ifelse(dsgmb2$ethnicity == "", NA, dsgmb2$ethnicity))

table(dsgmb2$ethnicity)

 Balanta     Fula     Jola Karonika Mandinka  Manjago    Other Sarahula 
       3       36       46        3       88        5       16        2 
   Serer   Wollof 
       5       25 
Display code
ethgmb <- dsgmb2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Gambia")+
  theme_bw()

ethgmb

Display code
## Pakistan:
eth <- as.data.frame(table(dspak$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dspak2 <- dspak[, c("Ethnicity")]

dspak2$ethnicity <- ifelse(dspak2$Ethnicity %in% l1, dspak2$Ethnicity, "Other")
dspak2$ethnicity <- ifelse(dspak2$ethnicity == "Serere", "Serer",
                 ifelse(dspak2$ethnicity == "", NA, dspak2$ethnicity))

table(dspak2$ethnicity)

      Bhatti    Christain       Hazara         Jopu        Other      Punjabi 
          22           28            7           20          124          118 
      Rajput      Sadozai         Shia         Sikh       Sindhi        Sunni 
           9           10           20           10           14           29 
        Syed Urduspeaking 
          34           60 
Display code
ethpak <- dspak2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Pakistan")+
  theme_bw()

ethpak

Display code
## Tanzania:
eth <- as.data.frame(table(dstza$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dstza2 <- dstza[, c("Ethnicity")]

dstza2$ethnicity <- ifelse(dstza2$Ethnicity %in% l1, dstza2$Ethnicity, "Other")
dstza2$ethnicity <- ifelse(dstza2$ethnicity == "Serere", "Serer",
                 ifelse(dstza2$ethnicity == "", NA, dstza2$ethnicity))

table(dstza2$ethnicity)

 Maasai  Mchaga   Mgogo Mlugulu Mluguru  Mngoni   Mpare Msambaa Msukuma    Muha 
    122       9       8      11      15       8       9      11      10       7 
 Mzigua   Other 
      7     120 
Display code
ethtza <- dstza2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Tanzania")+
  theme_bw()

ethtza

Display code
## Uganda:
eth <- as.data.frame(table(dsuga$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 7, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

dsuga2 <- dsuga[, c("Ethnicity")]

dsuga2$ethnicity <- ifelse(dsuga2$Ethnicity %in% l1, dsuga2$Ethnicity, "Other")
dsuga2$ethnicity <- ifelse(dsuga2$ethnicity == "Serere", "Serer",
                 ifelse(dsuga2$ethnicity == "", NA, dsuga2$ethnicity))

table(dsuga2$ethnicity)

     Adhola    Amudaama      Atesot       Bantu      Etesot      Itesot 
         18          11          61          11          56          31 
  Japadhola japhadollah Japhadollah     Munyole      Musoga     Nilotic 
         97          22          23           7           8          27 
      Other 
         76 
Display code
ethuga <- dsuga2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Uganda")+
  theme_bw()

ethuga

Section 2. Measures of interest

Variables of interest in the analysis:

  • Endorsement of Barrier Bound Leadership (BBL)

  • Endorsement of Barrier Crossing Leadership (BCL)

  • Ingroup fusion

  • Ingroup identification

  • Outgroup bonds (sum of outgroup fusion+identification)

  • Empathic concern

  • Perspective taking

  • Perceived history of discrimination

  • Positive contact with outgroup

  • Negative contact with outgroup

  • Imagistic items:

    • Event: Positive affect
    • Event: Negative affect
    • Event: Episodic recall
    • Event: Shared perception
    • Event: Reflection
    • Event: Transformative for individual
    • Event: Transformative for group
    • Imagistic event (sum of all above)
  • Perception of religion Freedom

  • Experience of religious freedom

Variable: Endorsement of Barrier Bound Leadership (BBL)

Display code
ds$bbl <- (ds$ENDBBL01+ ds$ENDBBL02+ ds$ENDBBL03)/3
ds$bcl <- (ds$ENDBCL01+ ds$ENDBCL02+ ds$ENDBCL03)/3

summary(ds$bbl)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       3       5       4       6       7      42 
Display code
ds %>% drop_na(bbl)%>%
ggplot(aes(x = bbl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(bbl, Country)%>%
ggplot(aes(x = bbl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(bbl)%>%
ggplot(aes(x = bbl,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Endorsement of BBL score", 
       y = "Frequency", 
       title = "Endorsement of BBL by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Endorsement of Barrier Crossing Leadership (BCL)

Display code
summary(ds$bcl)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      38 
Display code
ds %>% drop_na(bcl)%>%
ggplot(aes(x = bcl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(bcl, Country)%>%
ggplot(aes(x = bcl))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(bcl)%>%
ggplot(aes(x = bcl,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Endorsement of BCL score", 
       y = "Frequency", 
       title = "Endorsement of BCL by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Ingroup fusion

Display code
ds$igfusion <- (ds$IGF01+ds$IGF02+ds$IGF03)/3

summary(ds$igfusion)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      48 
Display code
ds %>% drop_na(igfusion)%>%
ggplot(aes(x = igfusion))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(igfusion, Country)%>%
ggplot(aes(x = igfusion))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(igfusion)%>%
ggplot(aes(x = igfusion,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup fusion score", 
       y = "Frequency", 
       title = "Ingroup fusion by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Ingroup identification

Display code
ds$IGI01 <- ifelse(ds$IGI01 %in% 1:7, ds$IGI01, NA)
ds$igidentification <- (ds$IGI01+ds$IGI02+ds$IGI03)/3

summary(ds$igidentification)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      48 
Display code
ds %>% drop_na(igidentification)%>%
ggplot(aes(x = igidentification))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(igidentification, Country)%>%
ggplot(aes(x = igidentification))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(igidentification)%>%
ggplot(aes(x = igidentification,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Ingroup identification score", 
       y = "Frequency", 
       title = "Ingroup identification by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Outgroup bonds

Display code
ds$ogbonds <- (ds$OGI01+ds$OGI02+ds$OGI03+
               ds$OGF01+ds$OGF02+ds$OGF03)/6

summary(ds$ogbonds)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       2       3       3       4       7     102 
Display code
ds %>% drop_na(ogbonds)%>%
ggplot(aes(x = ogbonds))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 3.00", 
                 xintercept = 3.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(ogbonds, Country)%>%
ggplot(aes(x = ogbonds))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(ogbonds)%>%
ggplot(aes(x = ogbonds,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 3.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Outgroup bonds score", 
       y = "Frequency", 
       title = "Outgroup bonds by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Empathic concern

Display code
ds$empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
                        ds$empathic_concern_03)/3

summary(ds$empathic_concern)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      33 
Display code
ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(empathic_concern, Country)%>%
ggplot(aes(x = empathic_concern))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Perspective taking

Display code
ds$perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
                          ds$perspective_taking_03+ds$perspective_taking_04)/4

summary(ds$perspective_taking)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      54 
Display code
ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 35)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(perspective_taking, Country)%>%
ggplot(aes(x = perspective_taking))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 25)+
 labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
 labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Perceived history of discrimination

Display code
summary(ds$history_discrimination)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       2       4       4       6       7      26 
Display code
ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 10)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(history_discrimination, Country)%>%
ggplot(aes(x = history_discrimination))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 10)+
 labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Positive Affect

Display code
summary(ds$Event_positive_affect)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       1       2       4       6       7      32 
Display code
ds %>% drop_na(Event_positive_affect)%>%
ggplot(aes(x = Event_positive_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 4.00", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_positive_affect, Country)%>%
ggplot(aes(x = Event_positive_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_positive_affect)%>%
ggplot(aes(x = Event_positive_affect,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 4.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Positive affect score", 
       y = "Frequency", 
       title = "Event: Positive affect by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Negative Affect

Display code
summary(ds$Event_negative_affect)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     2.0     6.0     4.8     7.0     7.0      26 
Display code
ds %>% drop_na(Event_negative_affect)%>%
ggplot(aes(x = Event_negative_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_negative_affect, Country)%>%
ggplot(aes(x = Event_negative_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_negative_affect)%>%
ggplot(aes(x = Event_negative_affect,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Negative affect score", 
       y = "Frequency", 
       title = "Event: Negative affect by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Episodic recall

Display code
summary(ds$Event_episodic_recall)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       6       6       6       7       7      42 
Display code
ds %>% drop_na(Event_episodic_recall)%>%
ggplot(aes(x = Event_episodic_recall))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_episodic_recall, Country)%>%
ggplot(aes(x = Event_episodic_recall))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_episodic_recall)%>%
ggplot(aes(x = Event_episodic_recall,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Episodic recall score", 
       y = "Frequency", 
       title = "Event: Episodic recall by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Shared perception

Display code
summary(ds$Event_shared_perception)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       6       7      32 
Display code
ds %>% drop_na(Event_shared_perception)%>%
ggplot(aes(x = Event_shared_perception))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_shared_perception, Country)%>%
ggplot(aes(x = Event_shared_perception))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_shared_perception)%>%
ggplot(aes(x = Event_shared_perception,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Shared perception score", 
       y = "Frequency", 
       title = "Event: Shared perception by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Reflection

Display code
ds$Event_reflection <- ds$Event_event_reflection
summary(ds$Event_reflection)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       4       6       5       6       7      40 
Display code
ds %>% drop_na(Event_reflection)%>%
ggplot(aes(x = Event_reflection))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_reflection, Country)%>%
ggplot(aes(x = Event_reflection))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_reflection)%>%
ggplot(aes(x = Event_reflection,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Reflection score", 
       y = "Frequency", 
       title = "Event: Reflection by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Transformative for individual

Display code
summary(ds$Event_transformative_indiv)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       4       6       5       6       7      49 
Display code
ds %>% drop_na(Event_transformative_indiv)%>%
ggplot(aes(x = Event_transformative_indiv))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_transformative_indiv, Country)%>%
ggplot(aes(x = Event_transformative_indiv))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_transformative_indiv)%>%
ggplot(aes(x = Event_transformative_indiv,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Transformative for individual score", 
       y = "Frequency", 
       title = "Event: Transformative for individual by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Event: Transformative for group

Display code
summary(ds$Event_transformative_group)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     4.0     5.0     4.9     6.0     7.0      28 
Display code
ds %>% drop_na(Event_transformative_group)%>%
ggplot(aes(x = Event_transformative_group))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_transformative_group, Country)%>%
ggplot(aes(x = Event_transformative_group))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_transformative_group)%>%
ggplot(aes(x = Event_transformative_group,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Transformative for group score", 
       y = "Frequency", 
       title = "Event: Transformative for group by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Variable: Imagistic event (sum of all subscales)

Display code
ds$Event_imagistic <- ds$Imagistic

summary(ds$Event_imagistic)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      7      32      36      35      40      49     102 
Display code
ds %>% drop_na(Event_imagistic)%>%
ggplot(aes(x = Event_imagistic))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
  geom_textvline(label = "Mean = 35.00", 
                 xintercept = 35.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic (full sample)")+
  theme_bw()

Display code
ds %>% drop_na(Event_imagistic, Country)%>%
ggplot(aes(x = Event_imagistic))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 15)+
 labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic by country")+
  facet_wrap(~Country)+
  theme_bw()

Display code
ds %>% drop_na(Event_imagistic)%>%
ggplot(aes(x = Event_imagistic,
           y = Country))+
  geom_boxplot(fill = "grey")+
  geom_textvline(label = " ", 
                 xintercept = 35.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+ 
  labs(x = "Event: Imagistic score", 
       y = "Frequency", 
       title = "Event: Imagistic by country")+
  #facet_wrap(~Country, nrow = 2)+
  theme_bw()

Section 3. Outcome: Social perception of Religious Freedom

Unconditional means model

Also called “varying intercept model with no predictors” (Gelman and Hill, 2016, Chapter 12). Allows intercepts to randomly vary across countries:

Display code
ds$SPRF <- ds$Sprf

## Varying intercept model with no predictors:
m00<- lmer(SPRF ~ 1 + (1 | Country), 
           data = ds)

summary(m00)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: SPRF ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 3068

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.726 -0.572  0.096  0.676  2.641 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0257   0.160   
 Residual             0.4502   0.671   
Number of obs: 1497, groups:  Country, 4

Fixed effects:
            Estimate Std. Error     df t value  Pr(>|t|)    
(Intercept)   5.4297     0.0822 3.0385      66 0.0000067 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:

Variance for Intercept = 0.03. This is the variance of the means across level 1 categories (countries). Residual = 0.44. Variance explained by level 1 residuals (everything that’s not in level 1).

Display code
tab_model(m00)
  SPRF
Predictors Estimates CI p
(Intercept) 5.43 5.27 – 5.59 <0.001
Random Effects
σ2 0.45
τ00 Country 0.03
ICC 0.05
N Country 4
Observations 1497
Marginal R2 / Conditional R2 0.000 / 0.054

We can see that ICC = 0.06. Lower ICC = low variance explained across groups. In this case, most of the variability is at individual-level (not group level). There is no significantly different patterns between countries.

Random intercept models

Also called “varying intercept model with individual-level predictors” (Gelman and Hill, 2016, Chapter 12).

Display code
## Create proper variable names:

ds$IG_Fusion <- (ds$IGF01+ds$IGF02+ds$IGF03)/3

ds$IG_Identification <- ds$igidentification

ds$OG_Bonds <- (ds$OGI01+ds$OGI02+ds$OGI03+
               ds$OGF01+ds$OGF02+ds$OGF03)/6

ds$Empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
                        ds$empathic_concern_03)/3

ds$Perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
                        ds$perspective_taking_03+ds$perspective_taking_03)/4

ds$Age <- ds$age

ds$Wealth_level <- as.factor(ifelse(ds$ses == "Lower middle", 1, 
                             ifelse(ds$ses == "Middle", 2,
                             ifelse(ds$ses == "Upper middle", 3,
                             ifelse(ds$ses == "Upper", 4, "")))))

ds$Female <- as.factor(ifelse(ds$gender == "Female", 1, 0))

ds$Married <- ds$married

ds$Perception_religious_freedom <- ds$SPRF

## Varying intercept models with individual-level predictors:
m01 <- lmer(Perception_religious_freedom~IG_Fusion+IG_Identification+OG_Bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m01)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Perception_religious_freedom ~ IG_Fusion + IG_Identification +  
    OG_Bonds + Empathic_concern + Perspective_taking + Age +  
    Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 2526

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-5.087 -0.593  0.033  0.677  4.070 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0192   0.139   
 Residual             0.3797   0.616   
Number of obs: 1315, groups:  Country, 4

Fixed effects:
                      Estimate  Std. Error          df t value
(Intercept)           2.791389    0.186736  126.440741   14.95
IG_Fusion             0.076776    0.022317 1301.703278    3.44
IG_Identification     0.124264    0.024302 1299.630921    5.11
OG_Bonds              0.017347    0.011905 1301.910330    1.46
Empathic_concern      0.038961    0.018456 1301.259994    2.11
Perspective_taking    0.174826    0.020026 1301.014809    8.73
Age                  -0.000787    0.001384 1301.732453   -0.57
Female1               0.039063    0.035127 1301.240598    1.11
MarriedOther         -0.015903    0.081643 1300.212099   -0.19
MarriedUnmarried      0.088790    0.040799 1300.996291    2.18
Wealth_level2         0.164497    0.044162 1201.020084    3.72
Wealth_level3         0.138814    0.063917 1244.155630    2.17
Wealth_level4        -0.035545    0.130529 1299.892380   -0.27
                               Pr(>|t|)    
(Intercept)        < 0.0000000000000002 ***
IG_Fusion                        0.0006 ***
IG_Identification            0.00000036 ***
OG_Bonds                         0.1453    
Empathic_concern                 0.0350 *  
Perspective_taking < 0.0000000000000002 ***
Age                              0.5694    
Female1                          0.2663    
MarriedOther                     0.8456    
MarriedUnmarried                 0.0297 *  
Wealth_level2                    0.0002 ***
Wealth_level3                    0.0301 *  
Wealth_level4                    0.7854    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Display code
tab_model(m01)
  Perception_religious_freedom
Predictors Estimates CI p
(Intercept) 2.79 2.43 – 3.16 <0.001
IG Fusion 0.08 0.03 – 0.12 0.001
IG Identification 0.12 0.08 – 0.17 <0.001
OG Bonds 0.02 -0.01 – 0.04 0.145
Empathic concern 0.04 0.00 – 0.08 0.035
Perspective taking 0.17 0.14 – 0.21 <0.001
Age -0.00 -0.00 – 0.00 0.569
Female [1] 0.04 -0.03 – 0.11 0.266
Married [Other] -0.02 -0.18 – 0.14 0.846
Married [Unmarried] 0.09 0.01 – 0.17 0.030
Wealth level [2] 0.16 0.08 – 0.25 <0.001
Wealth level [3] 0.14 0.01 – 0.26 0.030
Wealth level [4] -0.04 -0.29 – 0.22 0.785
Random Effects
σ2 0.38
τ00 Country 0.02
ICC 0.05
N Country 4
Observations 1315
Marginal R2 / Conditional R2 0.170 / 0.210

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code
## Change class of all models so we can use stargazer():
class(m00) <- "lmerMod"
class(m01) <- "lmerMod"

## Tabulated results:
stargazer(m00, m01,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")
Display code
htmltools::includeHTML("table1.html")
Dependent variable:
SPRFPerception_religious_freedom
(1)(2)
IG_Fusion0.077***
(0.022)
IG_Identification0.120***
(0.024)
OG_Bonds0.017
(0.012)
Empathic_concern0.039*
(0.018)
Perspective_taking0.170***
(0.020)
Age-0.001
(0.001)
Female10.039
(0.035)
MarriedOther-0.016
(0.082)
MarriedUnmarried0.089*
(0.041)
Wealth_level20.160***
(0.044)
Wealth_level30.140*
(0.064)
Wealth_level4-0.036
(0.130)
Constant5.400***2.800***
(0.082)(0.190)
Observations1,4971,315
Log Likelihood-1,534.000-1,263.000
Akaike Inf. Crit.3,074.0002,556.000
Bayesian Inf. Crit.3,090.0002,634.000
Note:*p<0.05; **p<0.01; ***p<0.001

Random intercept models: Imagistic predictors

Display code
## Varying intercept models with individual-level predictors:
m02 <- lmer(Perception_religious_freedom~Event_shared_perception+Event_episodic_recall+
              Event_reflection+Event_positive_affect+Event_negative_affect+
              Event_transformative_indiv+Event_transformative_group+
              Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m02)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Perception_religious_freedom ~ Event_shared_perception + Event_episodic_recall +  
    Event_reflection + Event_positive_affect + Event_negative_affect +  
    Event_transformative_indiv + Event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 2640

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.546 -0.647  0.025  0.683  3.091 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0359   0.189   
 Residual             0.3782   0.615   
Number of obs: 1369, groups:  Country, 4

Fixed effects:
                               Estimate   Std. Error           df t value
(Intercept)                   3.6878965    0.1683883   26.5844943   21.90
Event_shared_perception       0.0182817    0.0149981 1353.5600531    1.22
Event_episodic_recall         0.1463361    0.0180858 1352.8118832    8.09
Event_reflection              0.0071473    0.0147259 1343.1801914    0.49
Event_positive_affect         0.0127756    0.0098358 1353.0070788    1.30
Event_negative_affect         0.0343717    0.0094891 1353.3489940    3.62
Event_transformative_indiv    0.0448096    0.0152013 1353.6154566    2.95
Event_transformative_group    0.0369059    0.0139939 1352.2634171    2.64
Age                           0.0000562    0.0013581 1352.6211111    0.04
Female1                       0.0262406    0.0342730 1353.8540285    0.77
MarriedOther                 -0.0436935    0.0815275 1353.6452664   -0.54
MarriedUnmarried              0.0826492    0.0399443 1351.8882536    2.07
Wealth_level2                 0.1564298    0.0436726 1318.8752569    3.58
Wealth_level3                 0.0181246    0.0628422 1336.7050366    0.29
Wealth_level4                -0.1740043    0.1303406 1353.9648284   -1.33
                                       Pr(>|t|)    
(Intercept)                < 0.0000000000000002 ***
Event_shared_perception                 0.22308    
Event_episodic_recall        0.0000000000000013 ***
Event_reflection                        0.62751    
Event_positive_affect                   0.19421    
Event_negative_affect                   0.00030 ***
Event_transformative_indiv              0.00326 ** 
Event_transformative_group              0.00845 ** 
Age                                     0.96703    
Female1                                 0.44403    
MarriedOther                            0.59209    
MarriedUnmarried                        0.03873 *  
Wealth_level2                           0.00035 ***
Wealth_level3                           0.77307    
Wealth_level4                           0.18210    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Display code
tab_model(m02)
  Perception_religious_freedom
Predictors Estimates CI p
(Intercept) 3.69 3.36 – 4.02 <0.001
Event shared perception 0.02 -0.01 – 0.05 0.223
Event episodic recall 0.15 0.11 – 0.18 <0.001
Event reflection 0.01 -0.02 – 0.04 0.628
Event positive affect 0.01 -0.01 – 0.03 0.194
Event negative affect 0.03 0.02 – 0.05 <0.001
Event transformative
indiv
0.04 0.01 – 0.07 0.003
Event transformative
group
0.04 0.01 – 0.06 0.008
Age 0.00 -0.00 – 0.00 0.967
Female [1] 0.03 -0.04 – 0.09 0.444
Married [Other] -0.04 -0.20 – 0.12 0.592
Married [Unmarried] 0.08 0.00 – 0.16 0.039
Wealth level [2] 0.16 0.07 – 0.24 <0.001
Wealth level [3] 0.02 -0.11 – 0.14 0.773
Wealth level [4] -0.17 -0.43 – 0.08 0.182
Random Effects
σ2 0.38
τ00 Country 0.04
ICC 0.09
N Country 4
Observations 1369
Marginal R2 / Conditional R2 0.145 / 0.219

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code
## Change class of all models so we can use stargazer():
class(m00) <- "lmerMod"
class(m01) <- "lmerMod"
class(m02) <- "lmerMod"

## Tabulated results:
stargazer(m00, m01, m02,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")
Display code
htmltools::includeHTML("table1.html")
Dependent variable:
SPRFPerception_religious_freedom
(1)(2)(3)
IG_Fusion0.077***
(0.022)
IG_Identification0.120***
(0.024)
OG_Bonds0.017
(0.012)
Empathic_concern0.039*
(0.018)
Perspective_taking0.170***
(0.020)
Event_shared_perception0.018
(0.015)
Event_episodic_recall0.150***
(0.018)
Event_reflection0.007
(0.015)
Event_positive_affect0.013
(0.010)
Event_negative_affect0.034***
(0.009)
Event_transformative_indiv0.045**
(0.015)
Event_transformative_group0.037**
(0.014)
Age-0.0010.0001
(0.001)(0.001)
Female10.0390.026
(0.035)(0.034)
MarriedOther-0.016-0.044
(0.082)(0.082)
MarriedUnmarried0.089*0.083*
(0.041)(0.040)
Wealth_level20.160***0.160***
(0.044)(0.044)
Wealth_level30.140*0.018
(0.064)(0.063)
Wealth_level4-0.036-0.170
(0.130)(0.130)
Constant5.400***2.800***3.700***
(0.082)(0.190)(0.170)
Observations1,4971,3151,369
Log Likelihood-1,534.000-1,263.000-1,320.000
Akaike Inf. Crit.3,074.0002,556.0002,673.000
Bayesian Inf. Crit.3,090.0002,634.0002,762.000
Note:*p<0.05; **p<0.01; ***p<0.001

Histogram: Perception of Religious Freedom

Display code
summary(ds$SPRF)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      2       5       6       5       6       7      45 
Display code
ggplot(data = ds, 
       aes(x = SPRF)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.0", 
                 xintercept = 5.0, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "SPRF score", 
       title = "Social Perception of Religiour Freedom")+
  theme_bw()

Display code
ggplot(data = ds, 
       aes(x = SPRF)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "SPRF score", 
       title = "Social Perception of Religiour Freedom")+
  facet_wrap( ~ Country, nrow = 2) +
  theme_bw()

Display code
df01 <- ds %>% drop_na(SPRF)

tbl01 <- aggregate(df01$SPRF, 
                    by=list(df01$Country),
                    FUN=mean)
tbl01$Country <- tbl01$Group.1
tbl01$SPRF <- tbl01$x
tbl01 <- tbl01[, 3:4]
tbl01
   Country SPRF
1   Gambia  5.6
2 Pakistan  5.2
3 Tanzania  5.4
4   Uganda  5.6
Display code
ggplot(data = df01, 
       aes(x = SPRF, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 4.70", 
                 xintercept = 4.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "SPRF score", 
       title = "Social Perceptions of Relgious Freedom")+
  theme_bw()

Section 3. Outcome: Experience of Religious Freedom

Unconditional means model

Also called “varying intercept model with no predictors” (Gelman and Hill, 2016, Chapter 12). Allows intercepts to randomly vary across countries:

Display code
ds$Exp_religious_freedom <- ds$Exp_religious_freedom

## Varying intercept model with no predictors:
m10<- lmer(Exp_religious_freedom ~ 1 + (1 | Country), 
           data = ds)

summary(m10)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Exp_religious_freedom ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 5077

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.924 -0.515  0.369  0.705  1.210 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.129    0.36    
 Residual             1.743    1.32    
Number of obs: 1492, groups:  Country, 4

Fixed effects:
            Estimate Std. Error    df t value Pr(>|t|)    
(Intercept)    5.841      0.183 3.037    31.9 0.000062 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:

Variance for Intercept = 0.141. This is the variance of the means across level 1 categories (countries). Residual = 1.7. Variance explained by level 1 residuals (everything that’s not in level 1).

Display code
tab_model(m10)
  Exp_religious_freedom
Predictors Estimates CI p
(Intercept) 5.84 5.48 – 6.20 <0.001
Random Effects
σ2 1.74
τ00 Country 0.13
ICC 0.07
N Country 4
Observations 1492
Marginal R2 / Conditional R2 0.000 / 0.069

We can see that ICC = 0.08. Lower ICC = low variance explained across groups. In this case, most of the variability is at individual-level (not group level). There is no significantly different patterns between countries.

Random intercept models

Also called “varying intercept model with individual-level predictors” (Gelman and Hill, 2016, Chapter 12).

Display code
## Varying intercept models with individual-level predictors:
m11 <- lmer(Exp_religious_freedom~IG_Fusion+IG_Identification+OG_Bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m11)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Exp_religious_freedom ~ IG_Fusion + IG_Identification + OG_Bonds +  
    Empathic_concern + Perspective_taking + Age + Female + Married +  
    Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 4251

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.592 -0.524  0.233  0.696  2.405 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0852   0.292   
 Residual             1.4240   1.193   
Number of obs: 1316, groups:  Country, 4

Fixed effects:
                     Estimate Std. Error         df t value
(Intercept)           4.50051    0.36592   98.91990   12.30
IG_Fusion             0.04887    0.04321 1302.16380    1.13
IG_Identification     0.11225    0.04709 1300.48983    2.38
OG_Bonds             -0.20987    0.02305 1302.97689   -9.10
Empathic_concern      0.17408    0.03574 1302.85029    4.87
Perspective_taking   -0.00743    0.03884 1302.50507   -0.19
Age                   0.00260    0.00267 1302.57372    0.97
Female1               0.22350    0.06807 1302.67705    3.28
MarriedOther          0.22189    0.15809 1302.28377    1.40
MarriedUnmarried      0.04714    0.07906 1301.68725    0.60
Wealth_level2        -0.03748    0.08541 1235.25598   -0.44
Wealth_level3        -0.34497    0.12377 1264.82102   -2.79
Wealth_level4        -0.49288    0.25279 1302.08011   -1.95
                               Pr(>|t|)    
(Intercept)        < 0.0000000000000002 ***
IG_Fusion                        0.2583    
IG_Identification                0.0173 *  
OG_Bonds           < 0.0000000000000002 ***
Empathic_concern              0.0000013 ***
Perspective_taking               0.8484    
Age                              0.3303    
Female1                          0.0011 ** 
MarriedOther                     0.1607    
MarriedUnmarried                 0.5512    
Wealth_level2                    0.6609    
Wealth_level3                    0.0054 ** 
Wealth_level4                    0.0514 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Display code
tab_model(m11)
  Exp_religious_freedom
Predictors Estimates CI p
(Intercept) 4.50 3.78 – 5.22 <0.001
IG Fusion 0.05 -0.04 – 0.13 0.258
IG Identification 0.11 0.02 – 0.20 0.017
OG Bonds -0.21 -0.26 – -0.16 <0.001
Empathic concern 0.17 0.10 – 0.24 <0.001
Perspective taking -0.01 -0.08 – 0.07 0.848
Age 0.00 -0.00 – 0.01 0.330
Female [1] 0.22 0.09 – 0.36 0.001
Married [Other] 0.22 -0.09 – 0.53 0.161
Married [Unmarried] 0.05 -0.11 – 0.20 0.551
Wealth level [2] -0.04 -0.21 – 0.13 0.661
Wealth level [3] -0.34 -0.59 – -0.10 0.005
Wealth level [4] -0.49 -0.99 – 0.00 0.051
Random Effects
σ2 1.42
τ00 Country 0.09
ICC 0.06
N Country 4
Observations 1316
Marginal R2 / Conditional R2 0.122 / 0.172

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code
## Change class of all models so we can use stargazer():
class(m10) <- "lmerMod"
class(m11) <- "lmerMod"

## Tabulated results:
stargazer(m10, m11,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")
Display code
htmltools::includeHTML("table1.html")
Dependent variable:
Exp_religious_freedom
(1)(2)
IG_Fusion0.049
(0.043)
IG_Identification0.110*
(0.047)
OG_Bonds-0.210***
(0.023)
Empathic_concern0.170***
(0.036)
Perspective_taking-0.007
(0.039)
Age0.003
(0.003)
Female10.220**
(0.068)
MarriedOther0.220
(0.160)
MarriedUnmarried0.047
(0.079)
Wealth_level2-0.037
(0.085)
Wealth_level3-0.340**
(0.120)
Wealth_level4-0.490
(0.250)
Constant5.800***4.500***
(0.180)(0.370)
Observations1,4921,316
Log Likelihood-2,538.000-2,125.000
Akaike Inf. Crit.5,083.0004,281.000
Bayesian Inf. Crit.5,099.0004,359.000
Note:*p<0.05; **p<0.01; ***p<0.001

Random intercept models: Imagistic predictors

Display code
## Varying intercept models with individual-level predictors:
m12 <- lmer(Exp_religious_freedom~Event_shared_perception+Event_episodic_recall+
              Event_reflection+Event_positive_affect+Event_negative_affect+
              Event_transformative_indiv+Event_transformative_group+
              Age+Female+Married+Wealth_level+
              (1 | Country), 
            data = ds)

summary(m12)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Exp_religious_freedom ~ Event_shared_perception + Event_episodic_recall +  
    Event_reflection + Event_positive_affect + Event_negative_affect +  
    Event_transformative_indiv + Event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 4554

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.431 -0.583  0.280  0.722  1.891 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.126    0.355   
 Residual             1.571    1.253   
Number of obs: 1365, groups:  Country, 4

Fixed effects:
                             Estimate Std. Error         df t value
(Intercept)                   5.99486    0.33420   33.29321   17.94
Event_shared_perception       0.05464    0.03048 1349.92882    1.79
Event_episodic_recall         0.23886    0.03705 1349.09559    6.45
Event_reflection             -0.12720    0.03041 1329.62089   -4.18
Event_positive_affect        -0.09007    0.02018 1347.71554   -4.46
Event_negative_affect        -0.14057    0.01943 1349.59958   -7.23
Event_transformative_indiv    0.00229    0.03112 1349.11898    0.07
Event_transformative_group   -0.08194    0.02868 1348.51282   -2.86
Age                           0.00420    0.00276 1348.89116    1.52
Female1                       0.22681    0.06997 1349.95972    3.24
MarriedOther                  0.32713    0.16619 1349.86561    1.97
MarriedUnmarried              0.15721    0.08155 1348.03301    1.93
Wealth_level2                -0.14879    0.08869 1302.56527   -1.68
Wealth_level3                -0.45868    0.12821 1325.58354   -3.58
Wealth_level4                -0.68178    0.26560 1349.76810   -2.57
                                       Pr(>|t|)    
(Intercept)                < 0.0000000000000002 ***
Event_shared_perception                 0.07325 .  
Event_episodic_recall          0.00000000015806 ***
Event_reflection               0.00003074047694 ***
Event_positive_affect          0.00000871652381 ***
Event_negative_affect          0.00000000000078 ***
Event_transformative_indiv              0.94143    
Event_transformative_group              0.00434 ** 
Age                                     0.12849    
Female1                                 0.00122 ** 
MarriedOther                            0.04923 *  
MarriedUnmarried                        0.05409 .  
Wealth_level2                           0.09367 .  
Wealth_level3                           0.00036 ***
Wealth_level4                           0.01037 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Display code
tab_model(m12)
  Exp_religious_freedom
Predictors Estimates CI p
(Intercept) 5.99 5.34 – 6.65 <0.001
Event shared perception 0.05 -0.01 – 0.11 0.073
Event episodic recall 0.24 0.17 – 0.31 <0.001
Event reflection -0.13 -0.19 – -0.07 <0.001
Event positive affect -0.09 -0.13 – -0.05 <0.001
Event negative affect -0.14 -0.18 – -0.10 <0.001
Event transformative
indiv
0.00 -0.06 – 0.06 0.941
Event transformative
group
-0.08 -0.14 – -0.03 0.004
Age 0.00 -0.00 – 0.01 0.128
Female [1] 0.23 0.09 – 0.36 0.001
Married [Other] 0.33 0.00 – 0.65 0.049
Married [Unmarried] 0.16 -0.00 – 0.32 0.054
Wealth level [2] -0.15 -0.32 – 0.03 0.094
Wealth level [3] -0.46 -0.71 – -0.21 <0.001
Wealth level [4] -0.68 -1.20 – -0.16 0.010
Random Effects
σ2 1.57
τ00 Country 0.13
ICC 0.07
N Country 4
Observations 1365
Marginal R2 / Conditional R2 0.095 / 0.162

Here, marginal R sq is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code
## Change class of all models so we can use stargazer():
class(m10) <- "lmerMod"
class(m11) <- "lmerMod"
class(m12) <- "lmerMod"

## Tabulated results:
stargazer(m10, m11, m12,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")
Display code
htmltools::includeHTML("table1.html")
Dependent variable:
Exp_religious_freedom
(1)(2)(3)
IG_Fusion0.049
(0.043)
IG_Identification0.110*
(0.047)
OG_Bonds-0.210***
(0.023)
Empathic_concern0.170***
(0.036)
Perspective_taking-0.007
(0.039)
Event_shared_perception0.055
(0.030)
Event_episodic_recall0.240***
(0.037)
Event_reflection-0.130***
(0.030)
Event_positive_affect-0.090***
(0.020)
Event_negative_affect-0.140***
(0.019)
Event_transformative_indiv0.002
(0.031)
Event_transformative_group-0.082**
(0.029)
Age0.0030.004
(0.003)(0.003)
Female10.220**0.230**
(0.068)(0.070)
MarriedOther0.2200.330*
(0.160)(0.170)
MarriedUnmarried0.0470.160
(0.079)(0.082)
Wealth_level2-0.037-0.150
(0.085)(0.089)
Wealth_level3-0.340**-0.460***
(0.120)(0.130)
Wealth_level4-0.490-0.680*
(0.250)(0.270)
Constant5.800***4.500***6.000***
(0.180)(0.370)(0.330)
Observations1,4921,3161,365
Log Likelihood-2,538.000-2,125.000-2,277.000
Akaike Inf. Crit.5,083.0004,281.0004,588.000
Bayesian Inf. Crit.5,099.0004,359.0004,677.000
Note:*p<0.05; **p<0.01; ***p<0.001

Histogram: Experience of Religious Freedom

Display code
df01 <- ds%>% drop_na(Exp_religious_freedom)
summary(df01$Exp_religious_freedom)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     5.1     6.3     5.8     7.0     7.0 
Display code
ggplot(data = df01, 
       aes(x = Exp_religious_freedom)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.8", 
                 xintercept = 5.8, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "Exp_religious_freedom score", 
       title = "Experience of Relgious Freedom")+
  theme_bw()

Display code
ggplot(data = ds, 
       aes(x = Exp_religious_freedom)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "Exp_religious_freedom score", 
       title = "Experience of Religious Freedom")+
  facet_wrap( ~ Country, nrow = 2) +
  theme_bw()

Display code
df01 <- ds %>% drop_na(Exp_religious_freedom)

tbl01 <- aggregate(df01$Exp_religious_freedom, 
                    by=list(df01$Country),
                    FUN=mean)
tbl01$Country <- tbl01$Group.1
tbl01$Exp_religious_freedom <- tbl01$x
tbl01 <- tbl01[, 3:4]
tbl01
   Country Exp_religious_freedom
1   Gambia                   6.1
2 Pakistan                   5.4
3 Tanzania                   5.7
4   Uganda                   6.2
Display code
ggplot(data = df01, 
       aes(x = Exp_religious_freedom, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.8", 
                 xintercept = 5.8, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "Exp_religious_freedom score", 
       title = "Experience of Relgious Freedom")+
  theme_bw()

Section X: Coefficient Plots

Outcome: Social Perception of Religious Freedom

Display code
mp01 <- modelplot(m01,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Social perception of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp01

Display code
mp02 <- modelplot(m02,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Social perception of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp02

Display code
mp01a <- mp01 + theme(legend.position = "none")
ol2 <- mp01a / mp02
ol2

Outcome: Experience of Religious Freedom

Display code
mp11 <- modelplot(m11,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Experience of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp11

Display code
mp12 <- modelplot(m12,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Experience of religious freedom", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp12

Display code
mp11a <- mp11 + theme(legend.position = "none")
ol2a <- mp11a / mp12
ol2a